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A Systematic Description of Shocks in Gamma Ray Bursts 
II: Simulation 
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ABSTRACT 

In Paper I we presented a detailed formulation of the relativistic shocks and syn- 
chrotron emission in the context of Gamma-Ray Burst (GRB) physics. To see how 
well this model reproduces the observed characteristics of the GRBs and their after- 
glows, here we present the results of some simulations based on this model. They are 
meant to reproduce the prompt and afterglow emission in some intervals of time dur- 
ing a burst. We show that this goal is achieved for both short and long GRBs and 
their afterglows, at least for part of the parameter space. Moreover, these results are 
the evidence of the physical relevance of the two phenomenological models we have 
suggested in Paper I for the evolution of the active region, the synchrotron emitting re- 
gion in a shock. The dynamical active region model seems to reproduce the observed 
characteristics of prompt emissions and late afterglow better than the quasi-steady 
model which is more suitable for the onset of afterglows. Therefore these simulations 
confirm the arguments presented in Paper I about the behaviour of these models based 
on their physical properties. 

Key words: gamma-rays: bursts - Shockwaves. 



1 INTRODUCTION 

Gamma ray bursts show a variety of behaviours both in 
their early/prompt emissions in all energy bands and in 
their late emissions in lower energies - usually called the 
afterglow. There have been a large number of simulations to 
reproduce some of the aspects of GRB emission. Some of 
them are based on the phenomenological formulation of the 
evolu tion of a relativistic shock and its s ynchrotron emis- 
sion llPanaitescu et aZ.lll997l: iPanaitescrTfc Mcszaro^ 1 19981 : 
ILu ef a/.ll2006l : iDai fc Lull2002l : iMimica et a/., ,2008 ). Their 
main goal has been to show that the idea of shock origin 
for the GRBs is essentially correct. Others are based on 
heavy numerical magnetohydrodynamics simulations of rel- 
ativistic shocks and permit to reproduce some of the prop- 
erties of GRBs from fundamen tal principles (iLiang et al\ 



20061: iHededal fc Nordhmdl [200 5: Nishikaw a et aZ.I l2008l: 



Frederiksenll2008l:lNagataki et aZ.ll20071: ITchekhovskov ef al\ 



20081 ). But they present a small volume of the parameter 
space. None of these formulations and related simulations 
can be directly and systematically applied to real bursts 
data to estimate their parameters. 

With the simplified shock model explained in Paper I, 
we expect that the analytical approximation of the solutions 
of the kinematics and dynamics equations permit to simu- 
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late relativistic shocks and their synchrotron emission more 
consistently and realistically than previous simple formula- 
tions. Moreover, the relative simplicity of the formalism and 
the availability of analytical expressions for observables such 
as flux and lags should facilitate the extraction of the pa- 
rameters and the estimation of the evolution of important 
physical quantities. In its turn, this knowledge will help us 
to better understand the properties of the central engines of 
GRBs and their surroundings. 

A real GRB originates from one or more colliding 
shells. In Paper I we m entioned that simulations of the 
Fermi processes (jSpitkov skv 200^) show that the energy 
distribution of electrons and fields as well as the evo- 
lution of the synchrotron emission l|Nagataki et al\ l2007l : 
iTchekhovskov et al\ 1200^ 1 are complex and simple distri- 
butions such as power-law or even a power-law with ex- 
ponential cutoff are not good approximations for physical 
parameters. Nonetheless, within the precision of the avail- 
able data, they can be sufficiently close approximations in a 
short time interval. Therefore, we can model an entire burst 
by dividing it to separate regimes. However, we must first 
explore the parameter space and understand how various 
parameters affect the synchrotron emission that is our only 
observable. This task is not simple because the model has 
a large number of parameters. We leave the systematic ex- 
ploration of the parameter space to a future work. Here our 
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goal is to see if at least for some set of parameters this model 
can reproduce the observed features of the GRB emission. 

We consider a few set of parameters, determine various 
physical quantities explained in Paper I and their evolution, 
such as /3'(r'), a;^(r') (The definition of parameters are sum- 
marized in Table [5] and details can found in Paper I.), and 
observables such as light curves and spectrum. We remind 
that these simulations meant to correspond to a time inter- 
val in the life of a GRB and do not present the whole emis- 
sion. For determination of the light curves and spectrum we 
have used equation (60-Paper I) and we have neglected the 
terms proportional to G which are the contribution of high 
latitude emission and their Doppler shift. 



2 SIMULATIONS 

Figures [T] to [S] show the evolution of the quantities that de- 
termine the kinematics and dynamics of the shock and syn- 
chrotron emission. It is notable that they present a variety 
of behaviours despite the fact that in most cases parameters 
of the models are not too different from each others. Apriori 
some of these quantities can be extracted from observations 
(see Sec. 3 of Paper I). They give valuable information about 
the physical state of the ejecta and its surrounding. This in- 
formation can be used to constrain engine models. 

The first application of these plots is the comparison 
and verification of the present model against previous works. 
For instance, comparing plots of log (7' — 1) as a function of 



time in Fig. [3] with the simulations of Huang et al. (19 99) for 
relativistic and with simulations of iHuang fc Chend (|2003l ) 
for non-relativistic regime shows that although the general 
behaviour of this quantity in all simulations is the same, 
the slope of the evolution is not. According to the self- 
similar solution in the ultra-relativistic regime the slope of 

7'(t') 3/7 (iBlandford fc McKed Il977l : ISari et al\ 1 19981 : 

iHuang et aZ.II 19991 ). Here we find roughly the same value at 
the beginning of the emission, see specially models of after- 
glows in Fig. 13] where the relative Lorentz factor between two 
shells is highly relativistic. However, because in the models 
studied here the physical parameters such as the electric and 
magnetic fields are varying, the evolution of 7'(t') at the be- 
ginning of the shock can be much fiatter. On the other hand 
in non-relativistic regime, (7'(t') — 1) — > very quickly. This 
happens specially in the prompt shock where the relative 
Lorentz factor is small and approaches to zero when shells 
coalesce, see models 1 to 6 in Fig. [3] Therefore the trans- 
fer of kinematic energy to radiation is not really a simple 
power-law as usually assumed. This has important implica- 
tions for the formation of various features of the light curves 
(We discuss the light curves in details in the next section). 
There are very few systematic simulation of prompt emis- 
sion from first principles to compare with this results. There 
are howev er analytical estimations and simulations for the 
afterglow (|Huang et aZ.II 19991 : iHuang fc Cheng||2003l ). Their 
slope of a ~ —1 at very late times and somehow flatter 
at earlier earlier time is similar to models 9 and 10 which 
meant to simulate the late time afterglow. Nonetheless, the 
differences between these models as well as models 20 to 
22 which are also candidate models for late afterglow show 
that depending on the evolution of the active region and 
physical parameters, various kinematical and dynamical out- 



come can be expected from late afterglow. This is consistent 
with th e observation of very slow X-ray deca y of e.g. GRB 
050822 (|Blustin et fflZ.ll2005l : lGodet et al\2m^ with final de- 
cay slop e of g ~ —0.9 7 and ve ry fast decay o f e.g. GRB 
050730 (iHoUand et aZ.|[20 05: Pan dev et all l2006') and GRB 
060607a ifziaeepour et a/.i,2006b. . 120081) with final a < -3 

We notice that for models 17 to 19 which simulate the 
onset of prompt emission with an active region evolving 
according to quasi-steady state model, /?', 7', etc. have a 
pathological behaviour for one decade of time or so at the 
beginning of the simulation. We could not understand the 
reason because the same code has been used for the simula- 
tion of external shocks and they do not show any anormal 
behaviour. By changing the duration of the simulations and 
the same parameter sets we observed that the anormal in- 
terval varies, but the behaviour of normal section was stable 
and did not depend on the simulated duration. Therefore we 
conclude that the initial anormal behaviour is a purely nu- 
merical effect. Consequently, for these models we consider 
only the interval of time in which the simulation has good 
behaviour - from t — to ~ 0.1 sec. for models 17 and 18, and 
t — to ~ 0.03 sec. for model 19. 

Another physi cally important quantity is lo'^. 
In l|Sari et al\ 1 19981 ) its time evolution slope in radia- 
tive external shocks is estimated to be ~ —12/7 ~ —1.7. 
From plots in Fig. [5] we can conclude that in most cases 
the late time slope of is close to this value. However, 
we observe a more complex behaviour for this quantity. 
Notably, in prompt internal shocks it can increases at least 
during an interval of time. This is the case for models 1, 
17, 19, and 24. For most of other models the initial decline 
of is much slower than the predication of self-similar 
shocks. This means that the initial radiation in both 
internal and external shocks can be much harder than what 
is expected from self-similar shock models. This behaviour 
can be related to time variation of shocks micro-physics. 
We come back to this point in the next section. 

The column density of the active region n'Ar' is also an 
interesting quantity which apriori can be measured directly. 
Plots in Fig. [6] show that the intrinsic column density can 
evolve very differently from burst to burst, both in prompt 
internal and external shocks. Moreover its evolution depends 
on other physical parameters. In fact, from definition of this 
quantity it is clear that it directly depends on the evolution 
of density and the size of the active region see Figs. [7] and 
[8] A remarkable property is that during the sudden rise of 
the prompt emission such as in model 2, the column density 
simultaneously increases. The same in true during the build- 
up of active region in external shocks, models 13 and 15 to 
19. 

In the previ ous simulations l|Nakar fc PiranI 12004 
iKeshet et fflZ.|[200^ ) the effect of accumulation of matter has 
been added to the simulations as a factor called compression 
factor. Because in these simulations it is usually considered 
that synchrotron emission happens just in the shock front, 
the only column density seen by the radiation is considered 
to be the material in the line of sight and downstream of 
the shock front. Therefore we can not compare the column 
density of active zone and its evolution in this model with 
previous simulations. At present the extraction column den- 
sity of active region from observations is very difficult and 
depends on how precisely we can determine and remove the 
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contribution of foreground matter on the line of sight of 
GRBs. We have discussed these issues in Sec. 3 of Paper I. 



3 LIGHT CURVES 

The most prominent feature of the GRBs is the fast time 
variation in their emission reflected in the name of these 
transient sources. Any model of gamma ray bursts must re- 
produce their light curves and various observed features such 
as lags, Fast Rise Exponential Decay (FRED) behaviour of 
pulses, shallow slope regime in the early X-ray and some- 
times optical afterglow emission, chromatic breaks, and the 
slow decay of optical and longer wavelength emissions. 

To see if the model discussed in Paper I can reproduce 
the diversity of gamma-ray light curves and their afterglow 
in lower energy bands we calculate the time variation of the 
spectrum, equation (60-Paper I), assuming a power-law dis- 
tribution for the Lorentz factor of electrons. The parameter 
sets are the same as what we used for the calculation of the 
shell dynamics and emission characteristics summarized in 
Table (2] We should remind that multiple peaks observed in 
many GRBs are most probably the result of multiple shell 
collisions. Our aim here is to reproduce the characteristics of 
a single peak interpreted as the collision between two shells. 

We also remind that the parameters for the models with 
a dynamical active region have been chosen to simulate the 
prompt gamma-ray emission, except for models 7 to 10. For 
these models the parameter are chosen according to what 
we believe to correspond to the collision of the shells/ejecta 
with the low density slow material, with the wind around 
the central engine, or ISM. Such collisions are considered 
to be the origin of the late emissions that produces slowly 
varying emission at low energies. X-ray, optical, and radio. 
Although based on physical arguments the dynamical active 
region seems to be more suitable for the prompt emission, 
it can be also considered for the external shocks specially at 
late times. This is the aim of model 9 and 10. The chosen pa- 
rameters for the simulation of quasi steady models are most 
relevant to the afterglow, except models 17 to 19 for which 
they are similar to what we have chosen for the simulations 
of internal shocks with a dynamical active region. Similar to 
the previous case we want to see if the quasi steady model 
can reproduce the observed properties of the prompt emis- 
sion. The light curves obtained from these simulations are 
shown in Figs. [5] and 1101 

As mentioned in the introduction, in these simulations 
the high latitude emission - the terms proportional to G(r) 
defined in equation (61-Paper I) in the power spectrum, 
equation (60-Paper I) is neglected. This is a good approxi- 
mation for the prompt and early afterglow emissions. How- 
ever, at late times when the Lorentz factor is decreased, their 
contribution becomes important. The main effect of these 
retarded emission is a flattening of the lightcurve. Due to 
larger Doppler shift of higher latitudes with respect to the 
emission at zero latitude, it would be more significant in 
lower energies than higher ones. This effect is independent 
of the evolution of the active region and applies to all the 
models explained here. Therefore what ever the late slope of 
the simulations, we should assume that in full model they 
would be flatter. 



3.1 Light curve of simulated models with 
dynamical active region 

First we concentrate on the simulation of the models with a 
dynamical active region. Fig.|5]shows a variety of behaviours 
both in the time evolution of the light curves and in the dif- 
ference between the form of the light curve at different en- 
ergy bands. We should remind that as the value of Ar'(rQ) 
for a dynamical active region model can not be zero, we 
should assume that these simulations show only the evolu- 
tion of the collision after the formation of an active region 
that is used as the initial condition for these simulations. 
The plots in Fig. I^show that in some cases such as models 
3, 4, and 5 the parameter set does not permit an increase of 
the emission and it decays. By contrast in models 1, 2, and 
6 the initial emission evolves to a much harder and brighter 
emission. 

Model 1 which is hard and square-shape 
looks like some of the observed hard and very 



bright bursts such as G RB 060105 (IZiaeepour et al\ 
l2006al ). GRB 060813a (iMoretti et al\ l2006t ). GRB 



061007 jSchadv et al\ l2006b 



2007^ and the super burst 



GRB 080319b (|Racusin et aLM2008al lb[). The hardening is 
also reffected in the increasing 7^ and ui'm in this model 
see Figs. |4] and (5] Comparing the parameters of this model 
with other models, it is clear that the main difference is 
the evolution of the electric and magnetic fleld^. They 
both increase with radius, but the electric field rising 
slope is much sharper than magnetic field. Therefore it 
can maintain for relatively long time the acceleration 
of electrons. The growing magnetic field helps to have 
significant synchrotron emission but if the electric field was 
low, either all the synchrotron energy was emitted quickly 
- like model 2 - or the gross of emission was more gradual 
and made a triangle/Gaussian shape peak. The early decay 
of the emission in lowest energy bands is most probably 
due to numerical insufficiency. However, it is possible that 
the early fall of the optical and ultra-violet light curves are 
real and due to the hardening of the emission. 

In a real burst at some point the fields stop growing, 
the emission becomes softer and then breaks. Evidently, the 
evolution of fields according to a power-law with a constant 
index considered for each of the simulations here can not 
realistically reproduce the variation of parameters in a sin- 
gle simulation. Model 6, for instance, can be a prototype 
of long, few seconds, roughly constant emission following 
the rise of model 1. Model 2 shows a much faster rise sim- 
ilar to what happens in the short bursts. Fig. [3] shows that 
at maximum emission the shock becomes very quickly non- 
relativistic, i.e. 7' — > 1. One can predict that this should 
stop the growing of fields and in fact leads to their fast de- 
cay and thereby fast decay of the mission. This situation is 
well presented by model 3 in which the light curve decay in 
less than 0.1 sec. In this model the lag between high energy 
bands is ~ 20 msec. It should be possible to reduce the lag 
and find values consistent with what is observed for short 



^ When we talk about the evolution of the electric and magnetic 
fields, we actually mean the evolution of the fraction of kinetic 
energy of ejecta transferred to these fields. Fractions are scalar 
values in contrast to fields which are vectors and their variation 
is more complex. 
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bursts. Most probably this needs an exponential decay of 
the fields which is not included in the present work. Models 
4 and 5 are typical end of emission regime for long bursts in 
which a large lag, few hundreds of seconds, between energy 
bands appears. Another interesting aspect of this model is 
the very slow initial evolution of the optical and ultra-violet 
bands and relatively late rise of the optical bands. If the 
burst stays hard the peak of the emission in lower energy 
bands. X-ray and optical can have much larger lags. For in- 
stance, in models 4 and 5 soft X-ray and UV/optical bands 
are yet rising when the highest energy band in gamma-ray 
had decayed by more than one order of magnitude. Such 
large lag in low energy bands can be the origin of the peak 
obse rved in the X-ray - sp ecially in short hard bursts see 
e.g. (jZiaeepour et al H20q7al')- and optical light curve of some 
long bursts ( Oates et al\ 200^ ) . although an external shock 
origin can not be ruled out. On the other hand the steep 
decay slope of these models specially at softer bands can 
explain the very steep decline of the X-ray and optical emis- 
sions - usually called the tail emis sions- seen in the ma- 
jority of the Simft early afterglows (|Zhang et al\\200^ . In 
fact multi-band simultaneous observations o f late gamma- 
ray peaks in bursts such as GRB 0601 24 llHoUand et al\ 
2009 : iRomano et gZ.1 120081) , GRB 061121 (iPage et aZ.ll2006l . 



20071 1. and GRB 070724A l|Ziaeepour et aZ.ll2007al )liave al 



ready proved that this regime of X-ray light curves is in 
close relation with t he prompt emissio n (,Zhang et a/.ii2006 ; 
iKumar et oil l2006l : iLiang et al] 120071 ') consistent with the 
prediction of models presented here. In all these simulations 
the time duration of the emission (width of the peak) in 
higher energy bands are sh orter than lower energ ies, consis- 
tent with the observations (|Fenimore et aZ.ll 19951 ) 

Models 7 to 10 are simulated with parameters for af- 
terglows. From Fig. [8] it is clear that in model 7 the width 
of the active region is decreasing although it has a negative 
r. By contrast in model 8 which has a positive r the Ar' 
increases. Nonetheless in both cases the column density de- 
creases rapidly, see Fig. (6] These models have a relatively 
large Lorentz factor and hard emission in X-ray and their 
light curves are very flat or rising until hundreds of kilo 
seconds. Apriori this can explain the sh allow regime of the 
X-ray light curve s seen in many bursts (jZhang et aZ.|[200^ : 
iLiang et aZ.ll2007l ) if the shock stays only for a limited time 
of up to few thousands of seconds in this regime early af- 
ter the onset of the external shock. The late decay of light 
curves of these models is too fast and inconsistent with ob- 
servations. Therefore such a model can only be an example 
of early afterglows. We should also remember that the early 
emission in X-ray and optical band s is in part due to the 
tail emission of the prompt shock (|Willingale et aZ]|2007l ) 
and therefore a combination of the two component must be 
considered. We come back to this point in the next section. 

In model 9 and 10 all the parameters have values similar 
to the previous two models but their Lorentz factor is much 
smaller. In this case the shallow regime is declining, in con- 
trast to the previous models and the late decay has a slope of 
^ — 1.5 to —2.5 consistent with observations fZha ng et al\ 
I2OO6I ). The form of the light c urve also is gros sly consis- 
tent with previous simulations (|Sari et aZ.lll998l ) although 
the present models include more parameters. In fact if a 
regime similar to model 9 or 10 happens at late times > 10^ 
sec after trigger, in logarithmic scale the shallow part will 



be a very short and insignificant interval. Thus the observer 
sees a break from a previous regime, presumably a shallow 
regime similar to models 7 and 8, to a declining regime with 
a slope of ~ —1.5 to —2 that then breaks to a steeper slope 
of ~ —2.5 to ~ —3 at very late times. S uch a behaviour has 
been detected in many of Swift bursts (jO'Brien P.T. et al\ 
[2006!) . Our tests show that slopes somehow depend on p, 
larger p steeper the decli ne slope. T his is also consistent 
with previous simulations (|Sari et al\\l998,) . Nonetheless, as 
we explained in Paper I, flatter slopes can be obtained at 
low energies if there is an exponential cutoff in the electron 
distribution at high energies. This option is not included in 
present simulations and we leave its exploration to a future 
work. As mentioned before, the inclusion of high latitude 
emissions also make the decline flatter. 



3.2 Light curve of simulated models with 
quasi-steady active region 

Models 11 to 18 in Fig. [I0]are simulated with a quasi steady 
active region, the flrst six models as external shocks and the 
last three as internal shocks. The first noticeable feature is 
the similarity between the global trends of their light curves, 
in contrast to the variety of behaviours we have seen in mod- 
els with d ynamical active regi on. This is consistent with 



the Swift l|Gehrels et aZ.J |2 004) observations of X -ray and 
UV/optical late emissionfl (I Willingale et aLll2007h . Both X- 
ray and optical light curv es look like the lig ht curves of sim- 
ulated external shocks in ISari et al] (| 19981 ) . Their similarly 
to previous simulations is more visible than models 7 and 8 
with dynamical active region. 

In general these models show a large lag of few thou- 
sands of seconds between the time of peak emission at hard 
and soft X-ray bands, versus UV and opti cal bands. This is 
consistent with the previous calculations (|Sari et aZ.lll998l ). 
But the general aspects of the light curves are not sim- 
ilar to what is observed in real bursts. Therefore if the 
early emission was exclusively from external shocks these 
models were rules out. On t he other hand, it has been 
shown (|Willingale et a/.ll2007l ) that the X-ray light curves 
are composed of two overlapping components, one related 
to the prompt emission, and the other can be so called the 
afterglow. Moreover, the general form of the light curves 
in Fig. 1101 is similar and c onsistent with the fit performed 
in IWillingale et all (|2007l ) for the afterglow component. 
In this case the shallow regime observed in many bursts 
presents the interval in which the prompt compo nent is de- 
caying and the afterglow is gradually taking over (|loka et al\ 
I2OO6I : lo'Brien P. fc Willingalell200^ ). This can be also the 
reason for the hardening of the emission, because the decay- 
ing prompt emission can be at this time softer than growing 
emission from the external shock. As explained in the pre- 
vious section, another possibility for explaining the shallow 



^ We discussed that the observed emission in lower energies af- 
ter the prompt peak most probably consists of a combination of 
emission from the prompt collision and the afterglow defined as 
the emission from th e collision of ejecta/jet with the surround- 
ing material or ISM llWillingale et aZll2007l ). Therefore we refer 
to this emission as late emission rather than the afterglow which 
usually is considered to be the emission from external shocks 
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regime is the change in the evolution trend of the active re- 
gion from quasi-steady state to dynamical after the initial 
rise of emission in the previous regime. The slower decay 
rate of the low energy light curves in internal shocks and 
the large lag of the break in these bands with respect to X- 
ray break in the external shock emissions is consistent with 
the lack or very late break seen in the optical light curve s 
of the bursts detected by Swift /UVOT (Oatcs ei ai]|2008h . 
A shallower electron distribution at low energies and an ex- 
ponential breaJi at high ene rgy can also be re sponsible for 
shallow optical lightcurves (|Misra et al\\200^ . As for the 
decay slope after the break in models 11 to 13 it is about 
—2 and in model 14 about —2.5. Similar to dynamical mod- 
els for late external shock emission a shallower electron dis- 
tribution leads to a shallower break slope. In model 15 the 
break slope is very steep and inconsistent with observations. 
Therefore this model can simulate the rise of the afterglow 
and not its decay. Model 16 also has a steep slope but its 
Lorentz factor is not very large and high latitude emissions 
should significantly flatten the decay slope. 

In these simulations we have neglected the energy de- 
pendence of Ar, the width of the active region. The initial 
ejecta responsible for the internal shock is expected to be 
compact and dense. Thus, its emission in all bands come 
roughly from the same region and the assumption of an 
energy independent Ar is justified. Following the prompt 
collision(s) and scattering of the particles, one expects the 
formation of a low energy tail which nonetheless can have 
enough strong electric and magnetic fields to accelerate elec- 
trons and produce low energy synchrotron emission. To con- 
sider such a possibility we have also simulated the same 
models with the same parameters listed in Table [5] but with 
Aro Aro{Lij/ujm)~^^^ ■ The selected index of —1/2 for en- 
ergy dependence is arbitrary, and we leave a more precise 
estimation of this parameter to a future work in which we 
try to reconstruct real bursts according to the models pre- 
sented here. The results of these simulations are presented 
in Fig. 111! The general properties of the light curves are pre- 
served but lower energy bands have relatively higher fluxes 
and later breaks consistent with the theoretical explanation 
in Paper I. 

Models 17 to 19 are simulated with parameters relevant 
to internal shocks. They are generally hard and suitable for 
the rising regime of the prompt emission. In the dynamical 
active region model the initial Aro 7^ 0. Therefore somehow 
an initial active region must form before it evolves according 
to dynamical model. The quasi steady model is a suitable ap- 
proximation for this regime. The hardness of the simulations 
is consistent with the rapid hardening of the emission at the 
beginning of the collision. At high energies model 18 has a 
FRED shape light curve but the lag between gamma-ray en- 
ergy bands is ~ 2 sec. much longer than what is observed in 
long bursts, although very long faint bursts with supernova 
counterpart can accommodate such long lags. 

Finally Fig. [12] show the light curves for models 20 to 
24 in which the active region width evolves from an initial 
non-zero value at v'q according to a power-law, see equation 
(29-Paper I). The expression for Ar' evolution in this model 
apparently looks like quasi-steady state model, but in reality 
it is more similar to the dynamical model because in short 
intervals of time /3' and 7' on which the dynamical model de- 
pends, have a close to power-law evolution, see Figs. [2] and [S] 



Therefore, the evolution of Ar' is also a close to a power-law 
in a short time intervals. This fact is reflected in the form 
of the light curves of models 20 to 22 that present external 
shocks. They are very similar to models 7 to 10 generated by 
the dynamical model and their lightcurves have very shallow 
initial slopes and steep break. Their slopes after the break 
seem to be too steep. However, the contribution of the high 
latitude emissions as explained at the beginning of the Sec. 
[3] should flatten their decay slopes. Moreover, our tests and 
also the difference between the after-break slopes of models 
20 and 21 show that a flatter electron distribution reduces 
the steepness of the break. In particular, a power-law with 
exponential break distribution can make flatter late time 
breaks at lower energies and faster break at high energies. 
This feature has been also considered to exp lain the be- 
haviour of late X- ray omission of some bursts l|Misra et al\ 
2007i; Dempsey & Duflty 2007). Models 23 and 24 can be 
good candidate for the decay of emission during the prompt 
internal shock. 

In summary, the general behaviour of light curves in the 
simulations of active region models studied here are consis- 
tent with observations if we admit switching between models 
during the evolution of both internal and external shocks. 
Many of unexplained features of both prompt and afterglows 
are explained in these simulations such as: the energy depen- 
dence of the prompt emission peak width, the steep decay of 
the early afterglow and its connection to prompt emission, 
and the shallow regime which can be in part due to the way 
the active region evolves - in another word to the micro- 
physics of the shock - and in part due to the combination of 
tail emission from prompt shock and the onset of external 
shock. Anther notable conclusion from the simulations of ex- 
ternal shocks is that in all cases there is a chromatic break in 
the light curves due to the accumulation of in-falling matter 
into the shock. Therefore, changes of physical parameters 
and thereby the way the active region evolves, can com- 
pletely smear an achromatic break due to the kinematical 
jet break in a collimated ejecta. We observed a jet break in 
all the simulations although the model considered here has a 
spherical symmetry and therefore no kinematical jet break. 
Multiple jet breaks observed in a significant fraction of late 
emissions of the gamma-ray bursts also can be explained 
by switching between simple evolution models for the active 
region presented here as well as the change in the physical 
quantities such as the energy distribution of electrons. Ev- 
idently, it is up to magnetohydrodynamical simulations of 
relativistic and non-relativistic shocks and Fermi processes 
to verify the validity and conditions for the presence of the 
simple models presented here, the possibility of switching 
between models, different regimes, etc. 



4 SPECTRAL BEHAVIOUR 

In the GRB data analysis usually all the photons detected 
during a given interval of time, e.g. the duration of the ob- 
servation by the gamma-ray telescope - Tioo - is used for 
the determination of an averaged spectrum. Therefore for 
comparing the theoretical spectrum, equation (60-Paper I), 
with data we must integrate it with respect to r for a given 
radius/time interval and divide it by the duration length. 
However, here we do not consider a lower limit for the de- 
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tectable flux in the simulated models, and therefore Tioo is 
the same as the simulation time. For this reason we simply 
determine the total spectrum for the duration of the simu- 
lations without dividing it by time. The results are shown in 
Fig. [13] for dynamical and in Fig. [14] for quasi steady active 
region models explained in Table [2] 

4.1 Spectrum of simulated models with dynamical 
active region 

From variety of behaviours we have seen in the light curve 
of dynamical models we must expect the same level of vari- 
ation in the spectrum, and it is exactly what we find in Fig. 
1131 We should remind that the spectrum of the observed 
bursts are usually time averaged on the total duration of 
the burst, or belongs to the total duration of a peak. Thus, 
their behaviour can be very different from the spectrum in a 
shorter interval of time. Each of The simulations presented 
in this work correspond to a fraction of a burst, and there- 
fore we should not expect that it has a spectrum similar to 
time averaged observed spectrum of complete bursts. 

The spectra of models 1, 2, and 6 are much harder 
than observed time averaged spectra. However, considering 
their light curves, they correspond to a short time at the 
beginning of the burst, and therefor their impact on the 
integrated spectrum would be small. In f act some bright 
hard bursts such as G RB 060105 jZiaeepour et ali 
l2006al ). GRB 080319B jRacusin et ali l2008al l. and 
GRB 080706 l|Mang ano et al} l 2_008f) fl on g bursts) 



20051 : 



Fan fc Xul 
and GRB 



and GRB 051221 ( Parsons etaj^ 

l2006l ). GRB 060313 dPagani et al\ 120061 ') 
080123 (jUkwatta et ttZ.ll2008l 'l fshort bursts) seem to have a 
very flat or even rising spectrum at the onset of the burst. 

At highest energies, E > few times 10* eV, mod- 
els 3, 4, and 5 which present the end of a peak, see 
Fig. |§] have a power-law with an exponential cutoff spec- 
trum at i? > lO** eV and at lower energies the slope 
is positive. Th is is consistent with magnetohydrodynam- 
ics simulations (|Hededal fc Nordlund| l2005: Nishika wa et al\ 
[2OO8 ) and shock emission models |Medv edev 200^ . At 
E ^ 10* — 10^ eV where their spectra are roughly flat, 
with an error at the level or higher than what is shown 
in the plots (10% of the simulated values), they can be 
fit with a power-law consistent with observation of long 
Swift GRBs (Sakamoto et a/.i,2007. ). Although these models 
present only the end of a peak, the general behaviour of the 
broad band spectra is consistent the broad band spectrum of 
GRB 080916C which has been observed simultaneously by 



as Fermi llGoldstein & van der Horstll200S 


), Swift (iKennea 


20081: IStratta et 0Z.II2OO8I). Konus-Wind ( 


Golenetskii et al. 



20081 ). etc. This means that the effect of hard initial rise 



and soft tail emission after the decay of emission is small. 
Moreover, from similarity of these spectra with the spec- 
trum of the only burst for which we have broad spectrum 
up to GeV bands we can conclude that these models present 
typical prompt emissions. 

A common feature between all these spectra is their low 
luminosity at optical energies. This is consistent with the 
lack of very bright peak in optical frequencies when multiple- 
band observations of late g amma-ray peaks were possible. 
Examples are GRB 060526 (jCampana et aZ.ll2006l : rDai et al\ 



I2OO7I ) and GRB 061121 l|Page et aL|[2006l . |2007| ). In GRB 
060526 the X-ray emission increased about 100 times during 
a relatively soft but long peak observed also in gamma-ray. 
The optical emission during this peak increased only by a 
factor of < 2 in V band from the continuum. In GRB 061121 
also the X-ray increased about 100 times but V band peak 
was only ~ 5 times larger than the continuum and in B band 
only ~ 3 times. Most authors consider the null assumption 
of a single power-law spectrum at low energies and try to 
explain the low luminosity of bursts at low energies by dust 
absorption. In the case of GRB 061121 apriori a dust to gas 
model simi lar to the value for Milky Way can explain ob- 
servations (jPage et aLll200d ). However, we ignore the exact 
amount of dust in the host galaxy and we can not rule out a 
break at low en ergies consistent with magnetohydrodynam- 
ics simulations (jHededal fc Nordlundllioosi : iNishikawa et all 
2008). 

Models 7 to 10 simulate external shocks and have a 
much steeper spectrum at high energies. The spectrum of 
model 7 which can correspond to the rise of the external 
shock is hard with a break at few keV. The spectrum of 
model 8 and the soft X-ray spectra of models 9 and 10 can 
be fit by a single power-law with negative index. At very 
soft X and UV bands the spectrum index of models 9 and 

10 becomes steeper but at optical bands it becomes again 
flatter. This should correspond to the transition region from 
E > LJm to i? < ijjm where Um is the average minimum 
synchrotron characteristic energy in the observer frame. The 
spectrum indexes in X-ray band varies from ~ —0.4 in model 
8 to ~ -1.3 for model 9 and roughly flat for lO'^ < E < 10* 
eV for model 7. The large scatter of the X-ray spectrum 
index in these simu lations is consistent w ith the observations 
of the Swift -XRT ([Burrows et aLll2005l ). 

4.2 Spectrum of simulated models with 
quasi-steady active region 

Fig. [14] shows that the spectrum of external shock simula- 
tions according to this model is very similar to what have 
been obtained for the dynamical models. Except model 11 
which is hard and corresponds to the rising of the emission, 
others can be flt by a simple power-law or power-law with 
exponential cutoff at high and/or low energies. Like dynam- 
ical models, their spectrum index in X-ray band are mostly 
shallow with a > —1 and get even shallower in optical bands. 
But when the averaged minimum characteristic energy is in 
the plotted bands, as in model 12, there is a transient steep 
region similar to what we saw in dynamical models. Models 

11 and 15 are hard with a positive index at low energies. As 
these models can only be relevant for the rising of the ex- 
ternal shock emission which is usually overlapped with the 
tail emission of the prompt, it is very difficult to compare 
this results directly with observations. 

Fig. [TS] shows the spectrum of the same models but 
with an energy dependent Ar as explained in Sec l3.2l As 
expected, they are softer. Models 11 and 15 have negative 
slopes at low energies and a cutoff at high energies. Other 
models have the same form of the spectrum as in Fig. 1141 
but softer. 

Models 17 to 19 in Fig. [TJwhich present internal shocks 
with quasi-steady state active region are hard and has spec- 
tra similar to prompt models with dynamical active region. 



A Systematic Description of Shocks in Gamma Ray Bursts II: Simulation 7 



For various fundamental and observational reasons it 
is very difficult to compare these models with observations 
and with previous simulations or theoretical estimations. As 
we mentioned at the beginning of this section, these simu- 
lations does not correspond to a complete burst but just 
to part of a burst that can be simulated with a single and 
simple evolution model. Moreover, the emission after the 
prompt gamma-ray peak is most probably a combination of 
emission from t he prompt internal sho ck and the afterglow / 
external shocks (|Willingale et aZl2007l ) . It is not always easy 
to distinguish these components and in any case the separa- 
tion would be model dependent. Nonetheless, if we assume 
that the missing sectors do not change the spectrum signifi- 
cantly, we can compare these simulations with observations 
and previous theoretical estimations. For instance models 7, 
11, and 15 which present the early emission of the a fterglow 
with what is called fast cooling regime in dSari et a/.ii 19981 . 
the general form of the spectrum is similar. Notably, their 
slope at high energies is < — 1 when p = 2.5 for all models 
except model 16. This value i s con sistent with the value of 
—p/2 predicted in lSari et al\ (| 19981 ). In model 16 the index 
p = 3 and the value of slo pe at high energi es is ~ —1.5 
consistent with estimation of I Sari et al\ (|l998h . 

As for comparison with observations, it strongly de- 
pends on the time interval of observations. For instant, 
the available spectrum fit for the Swift bursts belong 
mainly to the first couple of thousands of seconds be- 
cause at later times the emission is faint and it is very 
difficult to determine the spectrum. However, as we men- 
tioned before the early X-ray and optical emission are 
most probably due to the tail emission from the rem- 
nant of the prompt shock and therefore they can not 
be compared with what is expected from external shock. 
This mixing can be the reason for the largo scatte r in 
the s pectral index of the X-ray spectrum (Willingalc et al. 
20071). For instance G RB 070724 A (Iziaeepour et al. 
2007al ). GRB 07072 IB JZiaeepour et al\ l2007bl '). GRB 



061004 jZiaeepour et al\ l2006d ) have an average spec- 
trum in dex of ^ —1.4 wh en others such as GRB 
05012 8 ((Cummings et g/.l 120051') . GRB 070208 jSato et al\ 
bOOTl ). and GRB 060804 (|Ziaeepour et aLll2006d ') have an 
index of ~ —2.3. 

On the other hand, one of the principle differences be- 
tween the prediction of these models and what is usually 
considered in the literature for the broad X-ray and opti- 
cal spectrum is a single power law with or without cutoff at 
high or low energies. The simulated models here show the 
possibility of variant slopes and presence of cutoffs even in 
the energy range accessible to one instrument. This can be 
the reason for the absence of an optical afterglow specially 
when the burst is hard. Nonetheless, the under-luminous 
optical emission can be also explained by the presence of 
dust i n the host galaxy or even in the circumburst environ- 
ment l|Schadv et a/.ll2007l ). One argument against the pres- 
ence of a large amount of dust in the host galaxy of GRBs 
is t he fact that most of them h ave relatively low metalic- 
ity (|Wolf fc Podsiadlowski|[2007l ). On the other hand, if the 
absorbers are limited to the star-forming clouds and their 
surroundings, then the local quantity of dust can be signif- 
icant although the total fraction of dust in the host galaxy 
can be low. Present observations and reduction techniques 
are not yet sufficient to clarify this issue. In conclusion, these 



simulations show that the spectrum of emission itself can be 
in part responsible for the lack of observation of optical and 
lower energy afterglows in the observed GRBs. 



5 SUMMARY 

We discussed the results of a few numerical simulations of 
GRBs based on the model explained in Paper I. We showed 
that the behaviour of these models is in general consistent 
with some regimes in the observed GRBs. It is difficult to 
verify the reality of some of the features seen in these sim- 
ulations. For instance, a long lag between the peak of the 
X-ray and optical afterglows is predicted by some of these 
simulations. We need long multi- wavelength observations to 
be able to see such a feature in real bursts, and usually such 
data do not exist. Nonetheless, there are exc eptions. Long 
follow -up of some bursts such as GRB 081029 (llmmler et ~al\ 
l2008h and GRB 081007 (jSoderberg et ffl;.ll2008l ) show a wide 
late peak in optical wavelengths undetected in the X-ray. It 
can be from the afterglow, although other possibility such as 
the appearance of a supernova can not be ruled out. These 
observations increase the confidence on this model and ana- 
lytical results presented here. Moreover, they encourage long 
term monitoring of GRBs to see if similar behaviour can be 
found in other bursts. 

To complete these simulations we need to put together 
the simulation of various regimes and make the light curves 
and spectra corresponding to the whole burst. This should 
help to better understand the contribution of internal and 
external shocks in the light curves and spectrum as well as 
features such as the shallow regime, breaks, etc. We leave 
these investigations for future works. 
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log |i'„,(r') 

-0.2 -0.1 




log fi'„(r') 
-0.4 -0.2 



log P'„(r-) 

-0.6 -0.4 -0.2 



log P'„(r') 
-0.15 -0.1 -0.05 



log »'jn 

-0.067 -0.066 -0.065 -0.064 







Pi 

n 
t3 



log fi'„(r') 

-0.075 -0.07 -0.065 



log P'„(r') 
-0.075 -0.07 -0.065 



1x10^-6x10^-4x10-^-2x10-^ 



log ?'„('■) 

-0.015 -0.01 -5x10-= 






log ?',,('•) 

-1.5x10-^ -10-^ 



log ?„{'•) 

-0.0635 -0.063 





log P'„(r') 
-0.064 -0.0635 -0.063 -0.0625 



log ^-„{f} 

-0.06 -0.04 -0.02 



log P'„(r') 
-0.04 -0.03 -0.02 -0.01 



log P'„(r') 
-0.15 -0.1 -0.05 



log P'„(r') 
-0.015 -0.01 -5x10-= 



log »',„{••) 

-0.0635 -0.063 
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t-t. sec 



105 
sec 




Figure 2. Evolution of /3j witli time/radius. Definitions of curves are the same as Fig. [T] The difference with /3g shows for each model 
the importance radiative correction in the evolution of the kinematics. 
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t-t g sec 

t-t g sec 

Figure 3. Evolution of the first order relative Lorentz factor with time/radius. In place of plotting 7J itself we plot log(7j — 1) to 
see more easily when the shock becomes non-relativistic. Definitions of curves are the same as Fig. ^ 
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Figure 4. Evolution of 7^ with time/radius. Definitions of curves are the same as Fig.[l] Note the hardening of emission with time for 
some of the models. 
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Figure 5. Evolution of ui!^ with timo/radius. Definitions of curves are the same as Fig.[T] These plots show that logui!^ and 7^ have a 
very similar evolution. 



A Systematic Description of Shocks in Gamma Ray Bursts II: Simulation 15 




21) 




22) 



23) 



24) 



1000 1C 
t-toSec 
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Figure 6. Evolution of the column density n'Ar' with time/radius. Definitions of curves are the same as Fig.[T] Depending on parameters 
and active region evolution model the column density showrs a variety of behavior. 
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t-t g sec 



Figure 7. Evolution of the active region density n' with time/radius. Definitions of curves are the same as Fig. ^ 
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Figure 8. Evolution of the width of the active region Ar' with time/radius. Definitions of curves are the same as Fig.^ 
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Table 1 . Definition of symbols 



Symbol 


Definition 




Aii c^ufintities witnout, except const3,nts, eltc QciineQ 




with respect to the rest frame of far but at the same 




redshift observer. 




rt.ll IJ U-dliLlLlCO WiLll d LJl illlC die tlcllllCtl WlLrli ItJOUCl^L LLJ 




the rest frame of slower shell. 


/5'(r') 


Relative fS of two colliding shell in the rest frame of 




slower shell. 




Characteristic synchrotron frequency for slowest elec- 




trons, Eq. (78-Paper I) 




Relative Lorentz factor of two shells in the rest frame 


of slower shell calculated with first order radiative cor- 




rection. 


n'Ar' 


intrinsic column density of the active region. 


P 


the index of electrons Lorentz factor distribution. 



Table 2. Parameter set of the simulated models. 





No. 


r'g {cm) 


Ar'„ 


P 


K 


7^ 


r 


<5 




as 






% (cm-3) 


n'K) 


^0 
































cm 


max 






1 


10" 


0.5 


2.5 





2 


-3 





10-3 


-0.5 


0.02 


-5 


5 X 10i2 


10" 


2 


200 




2 


5 X 10" 


0.1 


2.5 





2 


-3 





10-3 


-1 


0.3 


-3 


1013 


2 X 1013 


1.005 


200 




3 


5 X 10" 


0.2 


2.2 





2 


1 





2 X 10-3 


0.1 


0.2 





5 X 10i° 


2 X 10" 


1.005 


200 




4 


2 X 10" 


0.2 


3 





1.5 








10-3 


-0.3 


0.2 


-0.5 


5 X 10" 


1013 


1.1 


200 


D. 


5 


10" 


0.5 


2.7 





2 


-2 





10-4 


0.1 


0.2 


10-2 


5 X 10" 


1013 


1.5 


200 


A. 
R. 


6 


10" 


0.6 


2.1 





3 


-1 





10-3 


-0.8 


0.1 


-1 


5 Xl0i2 


5 X 1013 


3 


200 




7 


lO^b 


IQ-* 


2.5 


1 


50 


-3 





5 X 10-* 


0.5 


0.1 





1 


10'' 


5 


1 




8 


lOie 


10--* 


2.5 





50 


1 





5 X 10-3 


-0.01 


0.2 


-0.01 


50 


107 


5 


1 




9 


1016 


0.01 


2.5 


1 


10 


3 





5 X 10-3 


0.5 


0.1 





1 


IQS 


5 


1 




10 


1016 


0.01 


2.5 





10 


1 





5 X 10-3 


-0.01 


0.2 


-0.01 


10 


103 


5 


1 




11 


5 X 10i5 


0.05 


2.5 





200 





5 


10-3 


-0.5 


0.1 


-2 


10 


103 


5 


1 




12 


5 X 1016 


0.01 


2.5 


2 


20 





0.1 


10-3 


0.01 


0.5 


0.01 


10 


102 


2 


1 




13 


5 X 10i6 


0.01 


2.5 


2 


20 





1 


10-3 


0.01 


0.5 


0.01 


1 


10=2 


2 


1 


Q. 


14 


1015 


0.05 


2.5 





200 





5 


10-3 


-0.5 


0.1 


-1 


103 


10* 


10 


1 


S. 


15 


1016 


0.01 


2.5 


2 


200 





5 


10-3 


-1 


0.1 


-0.1 


1 


103 


5 


1 


A. 
R. - 


16 


1017 


5 xlO-3 


3 





10 





0.5 


10-3 


0.01 


0.2 


0.01 


10-1 


50 


2 


1 




17 


3 X 10^^ 


0.5 


2.5 





2 





4 


10-3 


0.01 


0.2 


-2 


lOl^ 


2 X 10^3 


1.5 


200 




18 


5 X 10" 


0.2 


2.5 





2 





4 


10-3 


0.2 


0.2 


2 


lOlO 


2 X 10" 


1.5 


200 




19 


10" 


0.5 


2.5 





2 





4 


10-3 


-0.5 


0.2 


-2 


1012 


2 X 1013 


1.5 


200 


Eq. 


20 


1016 


0.1 


2.5 





20 





0.5 


10-3 


0.01 


0.2 


0.01 


10 


103 


5 


1 


(29) 21 


1016 


0.1 


2.2 





10 





0.1 


10-3 


0.01 


0.2 


0.01 


10 


102 


2 


1 


in 


22 


5 X 10i6 


0.1 


2.5 





10 





1 


10-3 


0.05 


0.2 


0.05 


1 


102 


5 


1 


Pa- 


23 


3 X 10^^ 


0.5 


2.5 





2 





5 


10-3 


4 


0.2 


4 


1012 


2 X 10" 


1.5 


200 


per 


24 


5 X lOio 


0.5 


3 





2 





5 


10-3 


-0.5 


0.1 


-2 


1016 


2 X IQi** 


1.5 


200 



0Dynamical Active Region (D.A.R.), Quasi-Steady Active Region (Q.S.A.R.). 
Brhe aim for using large negative values for in the first and second parameter sets were simulating an exponential rise of the electric 



field. 
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1) 2) 3) 




10" 10= 10= 



t-t „ sec 

Figure 9. Light curves of models 1 to 10 in 10 energy bands from 2 eV to 2 MeV (models 1 to 6), 2 eV to 20 keV (models 7 and 8), 
and 2 eV to 2 keV (models 9 and 10), equally separated in logarithmic scale. The five highest energy bands (magenta/light grey), the 
five lowest energy bands (purple/dark grey). In each group the full line presents the highest energy band. Their active regions evolve 
according to dynamical model. Models 1 to 6 are meant to correspond to various regimes during prompt gamma-ray emission. Models 7 
and 8 presents the onset of afterglows, and models 9 and 10 the end of afterglows (See the text for detailed interpretation). 
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11) 12) 13) 




0.01 0.1 1 °'' ^ 0.01 0.1 1 



t-t„sec ' '"^^'^ t-t„sec 

Figure 10. Light curves of models 10 to 19 in 10 energy bands. For models 11 to 16 the energy bands are from 2 eV to 10 or 20 keV, 
equally separated in logarithmic scale. They are meant to correspond to onset of an afterglow. For models 17 to 19 the energy bands are 
from 2 eV to 2 MeV corresponding a prompt emission. Their active regions evolve according to quasi-steady state model. Descriptions 
of the curves are the same as Fig. |9] 
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11) 



14) 




13) 



2x10" 



15) 



Q. O r 





16) 




Figure 11. Light curves of models 11 to 16 in 10 energy bands with Ar oc iuj' juj'^ ^Z^. Descriptions of the energy bands and curves 
are the same as Fig. 1101 
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23) 




„ o 



24) 




100 1000 10« 
t-t „ sec 



o ^ 



0.01 



0.1 




Figure 12. Light curves of models 20 to model 24 in 10 energy bands, from 2 eV to 1 or 2 keV for models 20 to 22 and from 2 eV to 2 
MeV in models 23 and 24. Descriptions of the energy bands and curves are the same as Fig. 1101 Their active regions evolve according to 
equation (29-Paper I). Models 20 to 22 presents the end of of afterglow emission and models 23 and 24 the end of prompt emission. 
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Figure 13. Total broad band (2 cV - 2 McV)/(2 cV - 20 kcV) spectrum of the simulated models with dynamical active region. To show 
the level of ambiguity in the observational data due to systematic errors, we have added a random error with a normalized sigma of 5% 
for energy and 10% for spectrum amplitude (crosses). 



A Systematic Description of Shocks in Gamma Ray Bursts II: Simulation 23 



11) 12) 13) 




1000 10* 10= 10^ 10' ''° 1000 10" 10= 10« 

E (eV) ^ (^^' E (eV) 

Figure 14. Total broad band (2 eV - 20 keV) / (2 eV - 2 MeV) spectrum of the simulated models with quasi steady active region. 
Definitions of the curves are similar to Fig. 1131 
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11) 12) 13) 




Figure 15. Total broad band (2 eV - 20 keV) spectrum of the simulated models 11 to 16 with Ar oc (u'/'^m) Definitions of the 

curves are similar to Fig. 1131 
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10 100 1000 10* 10= 101= ^ 10 100 1000 10* 10= 10= 
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Figure 16. Total broad band (2 eV - 20 keV)/(2 eV - 2 MeV) spectrum of models 20 to 24. Definitions of the curves are similar to Fig. 

Ha 



